1. Summary

This report presents an exploratory analysis of social measures in Honduran communities participating in the Fish Forever (FF) program, focusing on trust, social cohesion, self-efficacy, and savings club participation. Due to key data limitations (absence of a control population, inconsistent temporal coverage at the community level, and low participation in savings clubs) the original objective of assessing program-driven change over time was revised. The analysis instead examines whether the duration of FF engagement is associated with current levels of the selected social measures. Regression results show no statistically significant relationships, with the exception of a marginal positive trend in perceived community ability to manage fisheries. Overall, the findings highlight critical limitations in the current dataset and point to the need for a more aligned and consistent sampling design if future survey efforts are intended to track change over time.

2. Original Objective

Request from the FF leadership:

“Find insights on social measures in Honduras (trust, social cohesion, self-efficacy). Specifically, assess if Fish Forever interventions are driving increases in these measures and if there are links to savings club participation”

3. Exploratory Anlaysis

The exploratory analysis is organized into two sections. The first, Available HHS Data, presents tables that summarize the household survey coverage in Honduras, disaggregated by province, municipality, and community across multiple years. The second, Selected HHS Questions, examines the response status and distribution of responses for a set of HHS questions aligned with the core social measures of interest: Trust, Social Cohesion, Self-Efficacy, and Savings Club Participation.

Available HHS Data

#Available Data
merged_hhs_hon <- merged_hhs %>% 
  filter(g1_country == "HND")


table_province <- merged_hhs_hon %>%
  group_by(year, merged_hhs_province) %>%
  summarise(n_surveys = n(), .groups = "drop")

table_municipality <- merged_hhs_hon %>%
  group_by(year, merged_hhs_province, merged_hhs_municipality) %>%
  summarise(n_surveys = n(), .groups = "drop")


table_community <- merged_hhs_hon %>%
  group_by(year, merged_hhs_province, merged_hhs_municipality, merged_hhs_community) %>%
  summarise(n_surveys = n(), .groups = "drop")


#Select relevant columns for analysis
merged_hhs_hon <- merged_hhs_hon %>%
  select(
    year,
    province = merged_hhs_province,
    municipality = merged_hhs_municipality,
    community = merged_hhs_community,
    g8_trust_community,
    g8_trust_community_emergency,
    g8_trust_community_neighbors,
    g8_trust_community_rebuilt,
    g8_trust_local_decision,
    g8_trust_regional_decision,
    g8_my_community_ability,
    g8_fishery_benefit_equal,
    g12_agreement_change_fishing_behavior,
    g12_agreement_individual_behavior,
    g7_savings_club_member,
    g7_emergency_savings_club 
  )


# Check missing data
# Define missing value patterns
missing_vals <- c(NA, 
                  "",
                  #"no_dependence", 
                  #"No dependance",
                  "Not Answered")

# Summarize missing and non-missing values per column
missing_summary <- map_dfr(names(merged_hhs_hon), function(col) {
  data <- merged_hhs_hon[[col]]
  missing_count <- sum(is.na(data) | data %in% missing_vals)
  total_count <- length(data)
  valid_count <- total_count - missing_count

  tibble(
    column = col,
    missing = missing_count,
    valid = valid_count,
    total = total_count,
    missing_pct = round(100 * missing_count / total_count, 1)
  )
})


# Province-level Table
table_province %>%
  arrange(merged_hhs_province, year) %>%
  rename(
    Province = merged_hhs_province,
    `Number of Surveys` = n_surveys,
    Year = year
  ) %>%
  arrange(Province) %>%
  select(Province, `Number of Surveys`, Year) %>%
  kable("html", caption = "Household Surveys by Province and Year") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Household Surveys by Province and Year
Province	Number of Surveys	Year
Atlántida	288	2023
Atlántida	288	2024
Colón	1160	2021
Colón	275	2023
Cortés	559	2019
Cortés	1	2024
Islas de la Bahía	227	2019



# Format table for display
table_municipality %>%
  arrange(merged_hhs_province, merged_hhs_municipality, year) %>%
  rename(
    Province = merged_hhs_province,
    Municipality = merged_hhs_municipality,
    `Number of Surveys` = n_surveys,
    Year = year
  ) %>%
  arrange(Province, Municipality) %>%
  select(Province, Municipality, `Number of Surveys`, Year) %>%
  kable("html", caption = "Household Surveys by Municipality and Year") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Household Surveys by Municipality and Year
Province	Municipality	Number of Surveys	Year
Atlántida	El Porvenir	288	2023
Atlántida	La Ceiba	288	2024
Colón	Balfate	3	2021
Colón	Balfate	275	2023
Colón	Iriona	185	2021
Colón	Limón	102	2021
Colón	Santa Fé	284	2021
Colón	Santa Rosa de Aguan	283	2021
Colón	Trujillo	303	2021
Cortés	Omoa	275	2019
Cortés	Puerto Cortés	284	2019
Cortés	Puerto Cortés	1	2024
Islas de la Bahía	Guanaja	227	2019



# Summarize years surveyed per community, including the actual years
community_years_surveyed <- table_community %>%
  group_by(merged_hhs_province, merged_hhs_municipality, merged_hhs_community) %>%
  summarise(
    `Years Surveyed` = paste0(
      n_distinct(year), 
      " (", paste(sort(unique(year)), collapse = ", "), ")"
    ),
    .groups = "drop"
  ) %>%
  arrange(merged_hhs_province, merged_hhs_municipality)

# Create the formatted table
community_years_surveyed %>%
  rename(
    Province = merged_hhs_province,
    Municipality = merged_hhs_municipality,
    Community = merged_hhs_community
  ) %>%
  kable("html", caption = "Years Surveyed per Community (Count and List of Years)") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Years Surveyed per Community (Count and List of Years)
Province	Municipality	Community	Years Surveyed
Atlántida	El Porvenir	Orotina	1 (2023)
Atlántida	El Porvenir	Porvenir Centro	1 (2023)
Atlántida	El Porvenir	Salado Barra	1 (2023)
Atlántida	La Ceiba	Corozal	1 (2024)
Atlántida	La Ceiba	La Ceiba Centro (Barra El Cangrejal)	1 (2024)
Atlántida	La Ceiba	La Ceiba Centro (Dantillo)	1 (2024)
Atlántida	La Ceiba	La Ceiba Centro (Miramar)	1 (2024)
Atlántida	La Ceiba	Sambo Creek	1 (2024)
Colón	Balfate	Balfate	1 (2023)
Colón	Balfate	Bambu	1 (2023)
Colón	Balfate	Lis Lis	1 (2023)
Colón	Balfate	Punta Gorda	1 (2021)
Colón	Balfate	Rio Coco	1 (2023)
Colón	Balfate	Rio Esteban	1 (2023)
Colón	Iriona	Cocalito	1 (2021)
Colón	Iriona	Cusuna	1 (2021)
Colón	Iriona	Iriona	1 (2021)
Colón	Iriona	Punta de Piedra	1 (2021)
Colón	Iriona	San José de la Punta	1 (2021)
Colón	Iriona	Sangre Laya	1 (2021)
Colón	Limón	Limón	1 (2021)
Colón	Santa Fé	Betulia	1 (2021)
Colón	Santa Fé	Guadalupe	1 (2021)
Colón	Santa Fé	Manatí	1 (2021)
Colón	Santa Fé	Plan Grande	1 (2021)
Colón	Santa Fé	Quinito	1 (2021)
Colón	Santa Fé	San Antonio	1 (2021)
Colón	Santa Fé	Santa Fe	1 (2021)
Colón	Santa Rosa de Aguan	Barra de Aguan	1 (2021)
Colón	Santa Rosa de Aguan	Col. Las Lomas	1 (2021)
Colón	Santa Rosa de Aguan	Dos Bocas	1 (2021)
Colón	Santa Rosa de Aguan	La Planada	1 (2021)
Colón	Santa Rosa de Aguan	Santa Rosa Aguan Centro	1 (2021)
Colón	Santa Rosa de Aguan	Vuelta Grande	1 (2021)
Colón	Trujillo	Barranco Blanco	1 (2021)
Colón	Trujillo	Capiro	1 (2021)
Colón	Trujillo	Cristales	1 (2021)
Colón	Trujillo	Jericó	1 (2021)
Colón	Trujillo	Marañones	1 (2021)
Colón	Trujillo	Moradel	1 (2021)
Colón	Trujillo	Puerto Castilla	1 (2021)
Colón	Trujillo	Rio Negro	1 (2021)
Colón	Trujillo	San Martín	1 (2021)
Colón	Trujillo	Silin	1 (2021)
Cortés	Omoa	Barra Motagua	1 (2019)
Cortés	Omoa	Las Flores	1 (2019)
Cortés	Omoa	Masca	1 (2019)
Cortés	Omoa	Milla 3	1 (2019)
Cortés	Omoa	Muchilena	1 (2019)
Cortés	Omoa	Omoa	1 (2019)
Cortés	Puerto Cortés	Bajamar	2 (2019, 2024)
Cortés	Puerto Cortés	Las Sabana	1 (2019)
Cortés	Puerto Cortés	Travesía	1 (2019)
Islas de la Bahía	Guanaja	El Cayo	1 (2019)
Islas de la Bahía	Guanaja	El Pelicano	1 (2019)
Islas de la Bahía	Guanaja	Mangrove Bight	1 (2019)
Islas de la Bahía	Guanaja	Savannah Bight	1 (2019)

Of the four provinces surveyed in Honduras, three have household survey data from more than one year. However, one of them (Cortés) includes only a single household surveyed in the second year (2024). When disaggregated to the community level, no community has data for more than one year (except Bajamar, which appears in both years but with only one household surveyed in 2024).

Selected HHS Questions

To assess the social outcomes intended by FF, we selected a subset of HHS questions that align with the core social measures of interest: Trust, Social Cohesion, Self-Efficacy, and Savings Club Participation. This section explores the response status (i.e., answered vs. missing) and the distribution of responses for each selected question. The analysis helps us understand both the availability and variability of data, serving as a foundation for the regression-based assessments that follow.

Trust

Community Trust: In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries.
Neighbors Trust: Generally speaking, most people in neighboring communities can be trusted.
Local Decision Makers: Local decision-makers/local authorities can be trusted to make decisions that benefit the community over their own interests.
Regional Decision Makers: Regional decision-makers/regional authorities can be trusted to make decisions that benefit the community over their own interests.

Social Cohesion

Community Support: My community will support me with funds and help in the case of emergency.
Post-Emergency Recovery: With the help of community members, the community can be rebuilt after emergencies occur.
Community Ability: My community has the ability to manage my fishery effectively to maximize food and profits.
Fishery Benefit Equality: Do you believe you benefit equally from the fishery as other members of the community?

Self-Efficacy

Change Behavior: I am willing to change my fishing behavior.
Individual Behavior: Through my individual fishing behavior, I can make a meaningful contribution to the sustainability of the fish catch.

Savings Club Participation

Are you or a member of your household a member of a savings club?
Does any household member have access to one of the following emergency funds? (Savings Club)

a. Trust

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_trust_community))
# 
# # Unique values
# unique(merged_hhs_hon$g8_trust_community)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_trust_community = case_when(
    str_detect(tolower(g8_trust_community), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_trust_community), "^disagree$|tidak setuju|b\\. tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_trust_community), "^neither$|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_trust_community), "^agree$|setuju|d\\. setuju") ~ "Agree",
    str_detect(tolower(g8_trust_community), "^strongly agree$|strongly_agree|sangat setuju|e\\. sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_trust_community), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))



################################

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_trust_community"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: c) Community Trust *",
    caption = "* c) In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )



# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_trust_community_neighbors))
# 
# # Unique values
# unique(merged_hhs_hon$g8_trust_community_neighbors)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_trust_community_neighbors = case_when(
    str_detect(tolower(g8_trust_community_neighbors), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju|sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_trust_community_neighbors), "^disagree$|b\\. tidak setuju|tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_trust_community_neighbors), "^neither$|neither agree nor disagree|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_trust_community_neighbors), "^agree$|d\\. setuju|setuju") ~ "Agree",
    str_detect(tolower(g8_trust_community_neighbors), "^strongly agree$|strongly_agree|e\\. sangat setuju|sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_trust_community_neighbors), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))



################################

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_trust_community_neighbors"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: d) Neighbors Trust *",
    caption = "* d) Generally speaking, most people in neighboring communities can be trusted",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_trust_local_decision))
# 
# # Unique values
# unique(merged_hhs_hon$g8_trust_local_decision)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_trust_local_decision = case_when(
    str_detect(tolower(g8_trust_local_decision), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju|sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_trust_local_decision), "^disagree$|b\\. tidak setuju|tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_trust_local_decision), "^neither$|neither agree nor disagree|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_trust_local_decision), "^agree$|setuju|d\\. setuju") ~ "Agree",
    str_detect(tolower(g8_trust_local_decision), "^strongly agree$|strongly_agree|sangat setuju|e\\. sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_trust_local_decision), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))




################################

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_trust_local_decision"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: a) Local Decision Makers *",
    caption = "* a) Local decision-makers/ local authorities can be trusted to make decisions that benefit the community over their own interests",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_trust_regional_decision))
# 
# # Unique values
# unique(merged_hhs_hon$g8_trust_regional_decision)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_trust_regional_decision = case_when(
    str_detect(tolower(g8_trust_regional_decision), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju|sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_trust_regional_decision), "^disagree$|b\\. tidak setuju|tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_trust_regional_decision), "^neither$|neither agree nor disagree|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_trust_regional_decision), "^agree$|setuju|d\\. setuju") ~ "Agree",
    str_detect(tolower(g8_trust_regional_decision), "^strongly agree$|strongly_agree|sangat setuju|e\\. sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_trust_regional_decision), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))




################################

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_trust_regional_decision"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: b) Regional Decision Makers *",
    caption = "* b) Regional decision-makers/ regional authorities can be trusted to make decisions that benefit the community over their own interests",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

b. Social Cohesion

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_trust_community_emergency))
# 
# # Unique values
# unique(merged_hhs_hon$g8_trust_community_emergency)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_trust_community_emergency = case_when(
    str_detect(tolower(g8_trust_community_emergency), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju|sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_trust_community_emergency), "^disagree$|b\\. tidak setuju|tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_trust_community_emergency), "^neither$|neither agree nor disagree|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_trust_community_emergency), "^agree$|setuju|d\\. setuju") ~ "Agree",
    str_detect(tolower(g8_trust_community_emergency), "^strongly agree$|strongly_agree|sangat setuju|e\\. sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_trust_community_emergency), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))




################################

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_trust_community_emergency"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(
    response_status == "Answered",
    response %in% response_levels,
    province %in% c("Atlántida", "Colón")
  ) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(c("Atlántida", "Colón")))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution for Atlántida and Colón",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: e) Community Support *",
    caption = "* e) My community will support me with funds and help in the case of emergency",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_trust_community_rebuilt))
# 
# # Unique values
# unique(merged_hhs_hon$g8_trust_community_rebuilt)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_trust_community_rebuilt = case_when(
    str_detect(tolower(g8_trust_community_rebuilt), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju|sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_trust_community_rebuilt), "^disagree$|b\\. tidak setuju|tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_trust_community_rebuilt), "^neither$|neither agree nor disagree|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_trust_community_rebuilt), "^agree$|setuju|d\\. setuju") ~ "Agree",
    str_detect(tolower(g8_trust_community_rebuilt), "^strongly agree$|strongly_agree|sangat setuju|e\\. sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_trust_community_rebuilt), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))




################################

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_trust_community_rebuilt"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(
    response_status == "Answered",
    response %in% response_levels,
    province %in% c("Atlántida", "Colón")
  ) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(c("Atlántida", "Colón")))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution for Atlántida and Colón",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: f) Post-Emergency Recovery *",
    caption = "* f) With the help of community members, the community can be rebuilt after emergencies occur",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )


# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_my_community_ability))
# 
# # Unique values
# unique(merged_hhs_hon$g8_my_community_ability)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_my_community_ability = case_when(
    str_detect(tolower(g8_my_community_ability), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju|sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_my_community_ability), "^disagree$|b\\. tidak setuju|tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_my_community_ability), "^neither$|neither agree nor disagree|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_my_community_ability), "^agree$|setuju|d\\. setuju") ~ "Agree",
    str_detect(tolower(g8_my_community_ability), "^strongly agree$|strongly_agree|sangat setuju|e\\. sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_my_community_ability), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))




################################

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_my_community_ability"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "How much do you agree with: “My community has the ability to manage my fishery effectively to maximize food and profits”?",
    #caption = "* c) In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_fishery_benefit_equal))
# 
# # Unique values
# unique(merged_hhs_hon$g8_fishery_benefit_equal)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_fishery_benefit_equal = case_when(
    str_detect(tolower(g8_fishery_benefit_equal), "^yes$|a\\. ya|^ya$") ~ "Yes",
    str_detect(tolower(g8_fishery_benefit_equal), "^no$|b\\. tidak|^tidak$") ~ "No",
    str_detect(tolower(g8_fishery_benefit_equal), "no_dependence") ~ "I don’t depend on or benefit from the fishery",
    str_detect(tolower(g8_fishery_benefit_equal), "not answered|^na$|tidak tahu|c\\. tidak tahu") ~ NA_character_,
    TRUE ~ NA_character_
  ))





################################

# Define missing values explicitly
missing_vals <- c(NA, 
                  "",
                  #"no_dependence", 
                  #"No dependance",
                  "Not Answered")

# Column to analyze
col <- "g8_fishery_benefit_equal"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Yes", "No", "I don’t depend on or benefit from the fishery"
)

color_mapping <- c(
"Yes" = "#6baed6",
    "No" = "#fb6a4a",
    "I don’t depend on or benefit from the fishery" = "#d9d9d9"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Do you believe you benefit equally from the fishery as other members of the community?",
    #caption = "* c) In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

c. Self Efficacy

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g12_agreement_change_fishing_behavior))
# 
# # Unique values
# unique(merged_hhs_hon$g12_agreement_change_fishing_behavior)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g12_agreement_change_fishing_behavior = case_when(
    str_detect(tolower(g12_agreement_change_fishing_behavior), "^strongly disagree$|^strongly_disagree$") ~ "Strongly disagree",
    str_detect(tolower(g12_agreement_change_fishing_behavior), "^disagree$") ~ "Disagree",
    str_detect(tolower(g12_agreement_change_fishing_behavior), "^neither agree nor disagree$|^neither$") ~ "Neither agree nor disagree",
    str_detect(tolower(g12_agreement_change_fishing_behavior), "^agree$") ~ "Agree",
    str_detect(tolower(g12_agreement_change_fishing_behavior), "^strongly agree$|^strongly_agree$") ~ "Strongly agree",
    str_detect(tolower(g12_agreement_change_fishing_behavior), "^$|not answered|^na$") |
      is.na(g12_agreement_change_fishing_behavior) ~ NA_character_,
    TRUE ~ NA_character_
  ))





################################

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g12_agreement_change_fishing_behavior"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: i) Change Behavior *",
    caption = "* i) I am willing to change my fishing behavior",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g12_agreement_individual_behavior))
# 
# # Unique values
# unique(merged_hhs_hon$g12_agreement_individual_behavior)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g12_agreement_individual_behavior = case_when(
    str_detect(tolower(g12_agreement_individual_behavior), "^strongly disagree$|^strongly_disagree$") ~ "Strongly disagree",
    str_detect(tolower(g12_agreement_individual_behavior), "^disagree$") ~ "Disagree",
    str_detect(tolower(g12_agreement_individual_behavior), "^neither agree nor disagree$|^neither$") ~ "Neither agree nor disagree",
    str_detect(tolower(g12_agreement_individual_behavior), "^agree$") ~ "Agree",
    str_detect(tolower(g12_agreement_individual_behavior), "^strongly agree$|^strongly_agree$") ~ "Strongly agree",
    str_detect(tolower(g12_agreement_individual_behavior), "^$|not answered|^na$") |
      is.na(g12_agreement_individual_behavior) ~ NA_character_,
    TRUE ~ NA_character_
  ))





################################


# Column to analyze
col <- "g12_agreement_individual_behavior"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: f) Individual Behavior *",
    caption = "* f) Through my individual fishing behavior, I can make a meaningful contribution to the sustainability of the fish catch",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

d. Savings Club

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g7_savings_club_member))
# 
# # Unique values
# unique(merged_hhs_hon$g7_savings_club_member)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g7_savings_club_member = case_when(
    str_detect(tolower(g7_savings_club_member), "^yes$") ~ "Yes",
    str_detect(tolower(g7_savings_club_member), "^no$") ~ "No",
    TRUE ~ NA_character_
  ))






################################


# Column to analyze
col <- "g7_savings_club_member"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Yes", "No"
)

color_mapping <- c(
"Yes" = "#6baed6",
    "No" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(
    response_status == "Answered",
    response %in% response_levels,
    province %in% c("Atlántida", "Colón")
  ) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(c("Atlántida", "Colón")))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution for Atlántida and Colón",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Are you or a member of your household a member of a savings club?",
    # caption = "* c) In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g7_emergency_savings_club))
# 
# # Unique values
# unique(merged_hhs_hon$g7_emergency_savings_club)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g7_emergency_savings_club = case_when(
    str_detect(tolower(g7_emergency_savings_club), "^yes$") ~ "Yes",
    str_detect(tolower(g7_emergency_savings_club), "^no$") ~ "No",
    str_detect(tolower(g7_emergency_savings_club), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))





################################


# Column to analyze
col <- "g7_emergency_savings_club"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Yes", "No"
)

color_mapping <- c(
"Yes" = "#6baed6",
    "No" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Does any household member have access to one of the following emergency funds? 6. Savings Club",
    #caption = "* c) In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

The plots show that several questions have significant missing data, indicating they were not included in all household surveys conducted in Honduras. Additionally, the plots reveal that Savings Club Participation is very low across surveyed communities, limiting the potential to explore meaningful relationships between savings club participation and social measures.

4. Limitations and Updated Objective

The original objective of this analysis, as outlined in Section 2, was to: “Find insights on social measures in Honduras (trust, social cohesion, self-efficacy). Specifically, assess if Fish Forever interventions are driving increases in these measures and if there are links to savings club participation.”

However, meeting this objective is not feasible due to three key limitations:

Lack of a Counterfactual (Control Population): The absence of non-intervention (control) communities prevents us from isolating the effects of FF on social measures. Any observed changes could be influenced by other factors such as economic shifts, governance changes, climate variability, or broader social dynamics; making attribution to FF alone unreliable.
Insufficient Temporal and Spatial Data Consistency: As shown in the tables above, no community has been surveyed in multiple years, making it impossible to assess changes over time at the community level. While a few provinces have been surveyed in more than one year, aggregating to that level would likely reflect geographic differences rather than true temporal change driven by program engagement.
Low Savings Club Participation: The plots above show that savings club participation is extremely limited across the surveyed sample. This lack of representation restricts the statistical power to explore meaningful relationships between participation and social outcomes.

Given these limitations, the revised objective of the analysis is to:

“Examine the relationship between the duration of Fish Forever’s engagement in communities and key social measures (trust, social cohesion, self-efficacy).”

This reframed objective focuses on cross-sectional associations between engagement length and current perceptions, rather than attempting to assess change or causal impact.

5. Analysis of FF Engagement Duration and Social Measures

program_start <- read_csv(here("data", "fp_level4_program_start.csv"), locale = locale(encoding = "ISO-8859-1"))

program_start <- program_start %>%
  mutate(level4_name = case_when(
    level4_name == "Dos Bocas / Brisas del Mar" ~ "Dos Bocas",
    level4_name == "Omoa casco" ~ "Omoa",
    level4_name == "Punta Piedra" ~ "Punta de Piedra",
    level4_name == "Iriona Viejo" ~ "Iriona",
    level4_name == "Estero / Las Flores" ~ "Las Flores",
    level4_name == "Santa Rosa Aguan Centro / La Planada" ~ "Santa Rosa Aguan Centro",
    level4_name == "Vuelta Grande/Fifi/Manguito" ~ "Vuelta Grande",
    level4_name == "Punta Gorda (JSG)" ~ "Punta Gorda",
    TRUE ~ level4_name
  ))

# Duplicate the row where level4_name == "Santa Rosa Aguan Centro"
duplicate_row <- program_start %>%
  filter(level4_name == "Santa Rosa Aguan Centro") %>%
  mutate(level4_name = "La Planada")

# Add it to the program_start dataframe
program_start <- bind_rows(program_start, duplicate_row)


# # Unique names from both datasets
# hhs_names <- unique(merged_hhs_hon$community)
# program_names <- unique(program_start$level4_name)
# # Find names in HHS not present in Program
# hhs_not_in_program <- setdiff(hhs_names, program_names)
# hhs_not_in_program

# Add program start to HHS data --> No program start data for Marañones
merged_hhs_hon <- merged_hhs_hon %>%
  left_join(
    program_start %>% select(level4_name, program_start),
    by = c("community" = "level4_name")
  )

merged_hhs_hon <- merged_hhs_hon %>% 
  filter(community != "Marañones")

In this analysis, we explore the potential impact of FF engagement duration on the three critical social measures: Trust, Social Cohesion, and Self-Efficacy. Using data collected through HHS across Honduran communities, we assess whether prolonged involvement in the Fish Forever program correlates with improvements in these measures. To enable this analysis, we incorporated data on when the program began in each community to calculate the duration of engagement. Linear regression models were then used to examine these relationships at the community level, accounting for variations in program duration across different sites.

a. Trust

# 1. Setup & Prepare the Data:

# Convert trust responses to numeric scale
likert_to_numeric <- function(x) {
  case_when(
    x == "Strongly disagree" ~ 1,
    x == "Disagree" ~ 2,
    x == "Neither agree nor disagree" ~ 3,
    x == "Agree" ~ 4,
    x == "Strongly agree" ~ 5,
    TRUE ~ NA_real_
  )
}

# Add numeric trust columns
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(
    trust_community_num = likert_to_numeric(g8_trust_community),
    trust_neighbors_num = likert_to_numeric(g8_trust_community_neighbors),
    trust_local_decision_num = likert_to_numeric(g8_trust_local_decision),
    trust_regional_decision_num = likert_to_numeric(g8_trust_regional_decision)
  )

# Convert dates and calculate duration of engagement
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(
    survey_date = make_date(year, 6, 30), # mid-year surveys
    program_start_date = dmy(program_start), # "18/4/2024" format
    ff_duration_years = as.numeric(difftime(survey_date, program_start_date, units = "days")) / 365,
    ff_duration_years = if_else(ff_duration_years < 0, 0, ff_duration_years) # negative durations as baseline (0)
  )

# 2. Aggregate Trust Scores at Community Level:

community_trust <- merged_hhs_hon %>%
  group_by(province, municipality, community) %>%
  summarise(
    mean_trust_community = mean(trust_community_num, na.rm = TRUE),
    mean_trust_neighbors = mean(trust_neighbors_num, na.rm = TRUE),
    mean_trust_local_decision = mean(trust_local_decision_num, na.rm = TRUE),
    mean_trust_regional_decision = mean(trust_regional_decision_num, na.rm = TRUE),
    avg_ff_duration_years = mean(ff_duration_years, na.rm = TRUE),
    n_surveys = n()
  ) %>%
  ungroup()

# 3. Linear Regression:

# Community Trust
lm_community_trust <- lm(mean_trust_community ~ avg_ff_duration_years, data = community_trust)

# Neighbors Trust
lm_neighbors_trust <- lm(mean_trust_neighbors ~ avg_ff_duration_years, data = community_trust)

# Local Decision-makers Trust
lm_local_decision_trust <- lm(mean_trust_local_decision ~ avg_ff_duration_years, data = community_trust)

# Regional Decision-makers Trust
lm_regional_decision_trust <- lm(mean_trust_regional_decision ~ avg_ff_duration_years, data = community_trust)

# Summary of Results
# summary(lm_community_trust)
# summary(lm_neighbors_trust)
# summary(lm_local_decision_trust)
# summary(lm_regional_decision_trust)

# 4. Visualization:

trust1 <- ggplot(community_trust, aes(x = avg_ff_duration_years, y = mean_trust_community)) +
  geom_point(size = 2, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "#1f78b4") +
  labs(
    title = "Community Trust",
    x = "Fish Forever Engagement Duration (Years)",
    y = "Avg. Community Trust"
  ) +
  theme_minimal()

trust2 <- ggplot(community_trust, aes(x = avg_ff_duration_years, y = mean_trust_neighbors)) +
  geom_point(size = 2, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "#1f78b4") +
  labs(
    title = "Neighbors Trust",
    x = "Fish Forever Engagement Duration (Years)",
    y = "Avg. Neighbors Trust"
  ) +
  theme_minimal()

trust3 <- ggplot(community_trust, aes(x = avg_ff_duration_years, y = mean_trust_local_decision)) +
  geom_point(size = 2, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "#1f78b4") +
  labs(
    title = "Local Decision Makers Trust",
    x = "Fish Forever Engagement Duration (Years)",
    y = "Avg. Local Decision-makers Trust"
  ) +
  theme_minimal()

trust4 <- ggplot(community_trust, aes(x = avg_ff_duration_years, y = mean_trust_regional_decision)) +
  geom_point(size = 2, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "#1f78b4") +
  labs(
    title = "Regional Decision Makers Trust",
    x = "Fish Forever Engagement Duration (Years)",
    y = "Avg. Regional Decision-makers Trust"
  ) +
  theme_minimal()


(trust1 + trust2) / (trust3 + trust4) +
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    title = "Relationship Between FF Engagement Duration and Trust",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20))
    )
  )

Interpretation of Results

Community Trust: There is no statistically significant relationship between FF engagement duration and community trust (p-value = 0.909). The coefficient is slightly negative (coefficient = -0.008), and the model has no explanatory power (Adjusted R² = -0.018).
Neighbors Trust: A small positive trend (coefficient = 0.062) exists between FF engagement duration and trust in neighboring communities, but it is not statistically significant (p-value = 0.415). The model explains very little of the variance (Adjusted R² = -0.006).
Local Decision-makers Trust: There is a mild positive association (coefficient = 0.085) between FF engagement duration and trust in local decision-makers, but the relationship is not statistically significant (p-value = 0.326). The model has almost no explanatory power (Adjusted R² = -0.0003).
Regional Decision-makers Trust: A similar mild positive trend (coefficient = 0.062) is observed for trust in regional authorities, though not statistically significant (p-value = 0.350). The adjusted R² is -0.002, indicating low explanatory value.

b. Social Cohesion

# 1. Likert to numeric conversion function
likert_to_numeric <- function(x) {
  case_when(
    x == "Strongly disagree" ~ 1,
    x == "Disagree" ~ 2,
    x == "Neither agree nor disagree" ~ 3,
    x == "Agree" ~ 4,
    x == "Strongly agree" ~ 5,
    x == "No" ~ 1,
    x == "Yes" ~ 5,
    TRUE ~ NA_real_
  )
}

# 2. Convert data and calculate FF duration
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(
    survey_date = make_date(year, 6, 30),
    program_start_date = dmy(program_start),
    ff_duration_years = as.numeric(difftime(survey_date, program_start_date, units = "days")) / 365,
    ff_duration_years = if_else(ff_duration_years < 0, 0, ff_duration_years)
  )

# 3. Social Cohesion Analysis
# Create numeric columns
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(
    community_support_num = likert_to_numeric(g8_trust_community_emergency),
    post_emergency_recovery_num = likert_to_numeric(g8_trust_community_rebuilt),
    community_ability_num = likert_to_numeric(g8_my_community_ability),
    fishery_benefit_equal_num = likert_to_numeric(g8_fishery_benefit_equal)
  )

# Aggregate at community level
social_cohesion <- merged_hhs_hon %>%
  group_by(province, municipality, community) %>%
  summarise(
    mean_community_support = mean(community_support_num, na.rm = TRUE),
    mean_post_emergency_recovery = mean(post_emergency_recovery_num, na.rm = TRUE),
    mean_community_ability = mean(community_ability_num, na.rm = TRUE),
    mean_fishery_benefit_equal = mean(fishery_benefit_equal_num, na.rm = TRUE),
    avg_ff_duration_years = mean(ff_duration_years, na.rm = TRUE),
    n_surveys = n()
  ) %>% ungroup()

# Linear regression analyses
lm_community_support <- lm(mean_community_support ~ avg_ff_duration_years, data = social_cohesion)
lm_post_emergency_recovery <- lm(mean_post_emergency_recovery ~ avg_ff_duration_years, data = social_cohesion)
lm_community_ability <- lm(mean_community_ability ~ avg_ff_duration_years, data = social_cohesion)
lm_fishery_benefit_equal <- lm(mean_fishery_benefit_equal ~ avg_ff_duration_years, data = social_cohesion)

# Summarize regression results
# summary(lm_community_support)
# summary(lm_post_emergency_recovery)
# summary(lm_community_ability)
# summary(lm_fishery_benefit_equal)

# Visualize results
plot_cohesion <- function(df, y_var, y_label, title) {
  ggplot(df, aes(x = avg_ff_duration_years, y = .data[[y_var]])) +
    geom_point(size = 2, alpha = 0.7) +
    geom_smooth(method = "lm", se = TRUE, color = "#1f78b4") +
    labs(
      title = title,
      x = "Fish Forever Engagement Duration (Years)",
      y = y_label
    ) +
    theme_minimal()
}

plot_cohesion1 <- plot_cohesion(social_cohesion, "mean_community_support", "Avg. Community Support", "Community Support")
plot_cohesion2 <- plot_cohesion(social_cohesion, "mean_post_emergency_recovery", "Avg. Post-Emergency Recovery", "Post-Emergency Recovery")
plot_cohesion3 <- plot_cohesion(social_cohesion, "mean_community_ability", "Avg. Community Fishery Management Ability", "Community Ability")
plot_cohesion4 <- plot_cohesion(social_cohesion, "mean_fishery_benefit_equal", "Avg. Equal Benefit from Fishery", "Fishery Equal Benefit")


(plot_cohesion1 + plot_cohesion2) / (plot_cohesion3 + plot_cohesion4) +
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    title = "Relationship Between FF Engagement Duration and Social Cohesion",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20))
    )
  )

Interpretation of Results

Community Support: There is no statistically significant relationship between FF engagement duration and perceptions of community support (p-value = 0.723). The coefficient is small (coefficient = 0.034), and the model explains virtually no variation (Adjusted R² = -0.072).
Post-Emergency Recovery: Likewise, there is no significant association between FF duration and perceptions of post-emergency recovery (p-value = 0.789). The coefficient is slightly negative (coefficient = -0.017), and the adjusted R² is -0.077.
Community Ability (marginal result): There is a slight positive trend between FF engagement duration and perceived ability of the community to manage fisheries (coefficient = 0.125), with marginal significance (p-value = 0.095). The model has modest explanatory power (Adjusted R² = 0.033).
Fishery Benefit Equality: No relationship is found between FF engagement duration and perceived equality in fishery benefits (p-value = 0.932). The coefficient is nearly zero (coefficient = -0.007), and the adjusted R² is -0.018.

c. Self Efficacy

# 4. Self-Efficacy Analysis
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(
    change_behavior_num = likert_to_numeric(g12_agreement_change_fishing_behavior),
    individual_behavior_num = likert_to_numeric(g12_agreement_individual_behavior)
  )

# Aggregate at community level
self_efficacy <- merged_hhs_hon %>%
  group_by(province, municipality, community) %>%
  summarise(
    mean_change_behavior = mean(change_behavior_num, na.rm = TRUE),
    mean_individual_behavior = mean(individual_behavior_num, na.rm = TRUE),
    avg_ff_duration_years = mean(ff_duration_years, na.rm = TRUE),
    n_surveys = n()
  ) %>% ungroup()

# Linear regression analyses
lm_change_behavior <- lm(mean_change_behavior ~ avg_ff_duration_years, data = self_efficacy)
lm_individual_behavior <- lm(mean_individual_behavior ~ avg_ff_duration_years, data = self_efficacy)

# Summarize regression results
# summary(lm_change_behavior)
# summary(lm_individual_behavior)

# Visualize results
plot_efficacy <- function(df, y_var, y_label, title) {
  ggplot(df, aes(x = avg_ff_duration_years, y = .data[[y_var]])) +
    geom_point(size = 2, alpha = 0.7) +
    geom_smooth(method = "lm", se = TRUE, color = "#1f78b4") +
    labs(
      title = title,
      x = "Fish Forever Engagement Duration (Years)",
      y = y_label
    ) +
    theme_minimal()
}

plot_efficacy1<- plot_efficacy(self_efficacy, "mean_change_behavior", "Avg. Willingness to Change Behavior", "Willingness to Change Behavior")
plot_efficacy2 <- plot_efficacy(self_efficacy, "mean_individual_behavior", "Avg. Belief in Individual Impact", "Belief in Individual Impact")


(plot_efficacy1 / plot_efficacy2)  +
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    title = "Relationship Between FF Engagement Duration and Self Efficacy",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20))
    )
  )

Interpretation of Results

Change Behavior: There is no significant relationship between FF engagement duration and willingness to change fishing behavior (p-value = 0.739). The coefficient is small (coefficient = 0.021), and the model has no explanatory value (Adjusted R² = -0.016).
Individual Behavior Impact: A small positive trend exists between FF engagement and belief in one’s individual impact on sustainability, but the result is not statistically significant (p-value = 0.220). The coefficient is small (coefficient = 0.094), with low explanatory power (Adjusted R² = 0.010).

6. Conclusions

Due to key limitations in the available data (absence of a control population, a lack of temporally consistent data at the community level, and very limited participation in savings clubs) we revised the original objective of this analysis. Rather than assessing whether FF interventions have caused improvements in social measures over time or are linked to savings club participation, we focused instead on examining whether the duration of FF engagement in a community is associated with current levels of trust, social cohesion, and self-efficacy.

Overall, the regression analyses indicate limited evidence of such a relationship. Trust metrics showed no statistically significant association with the length of FF engagement. Likewise, indicators of social cohesion and self-efficacy did not exhibit significant correlations with program duration. One exception was a weak significant positive relationship between FF engagement duration and perceived community ability to manage fisheries, which may suggest some marginal incremental gains over time.

Ultimately, this report highlights the limitations of the current dataset more than it reveals actionable findings. The fragmented and inconsistent structure of the data poses serious constraints for assessing change over time or attributing impact. With these limitations comes a clear recommendation: if the goal of future survey efforts is to track changes in social measures, FF will need to adopt an aligned and consistent sampling design. Notably, I understand that FF is already aware of these challenges and is considering improvements in this area, which will be essential for generating more meaningful and evaluable data in future assessments.

Social Measures in Honduras: Data Insights and Limitations

Mariano Viz

2025-07-08

1. Summary

2. Original Objective

3. Exploratory Anlaysis

Available HHS Data

Selected HHS Questions

a. Trust

c. Self Efficacy

d. Savings Club

4. Limitations and Updated Objective

6. Conclusions